#!/usr/bin/python

import sys
import re
from Bio.Blast import NCBIXML
#from Bio import SeqIO

xmlfile = sys.argv[1]
#seqfile = sys.argv[2]

result_handle = open(xmlfile)
#sf = SeqIO.read(seqfile, "fasta")

m1 = re.compile('\w+_\w+_\d+_(\d+)-(\d+).*')
m2 = re.compile('^PREDICTED.*', re.IGNORECASE)
m3 = re.compile('^hypothetical.*', re.IGNORECASE)

pattern = m1.match(xmlfile)

igs_start = pattern.group(1)
igs_stop = pattern.group(2)

blast_record = NCBIXML.read(result_handle)

desc_list = []
def_list = []
fstart_list = []
fstop_list = []
query_len_list = []
hit_len_list = []
#seq_list = []
hsp_list = []
hsp_starts = []
hsp_stops = []
top_hsp_starts = []
top_hsp_stops = []

for hsps in blast_record.descriptions:
    hsp_list.append(hsps.num_alignments)

#print hsp_list[0]

for alignment in blast_record.alignments:
    for hsp in alignment.hsps:
        hsp_starts.append(hsp.query_start)
        hsp_stops.append(hsp.query_end)
        #fstart_list.append(feat_start)
        #fstop_list.append(feat_stop)
        def_list.append(alignment.hit_def)
        #aa_length = (feat_stop - feat_start + 1) / 3
        #query_len_list.append(aa_length)
        #seq = sf.seq[feat_start:feat_stop]
        #seq_list.append(seq)
        hit_len_list.append(alignment.length)

#print seq_list[0], len(seq_list[0])
#print seq_list[0].reverse_complement().translate(), len(seq_list[0].translate())

num = int(hsp_list[0])

i = 0

while i < num:
    top_hsp_starts.append(hsp_starts[i])
    top_hsp_stops.append(hsp_stops[i])
    #print hsp_starts[i], hsp_stops[i]
    i += 1

top_hsp_starts.sort()
top_hsp_stops.sort()
smallest_start = top_hsp_starts[0]
largest_stop = top_hsp_stops[num - 1]

feat_start = int(igs_start) + int(smallest_start) - 1
feat_stop = int(igs_start) + int(largest_stop)
aa_length = (feat_stop - feat_start + 1) / 3

if len(def_list) != 0:
    if m2.match(def_list[0]) or m3.match(def_list[0]):
        pass
    else:
        if int(aa_length) < float(hit_len_list[0]) * 0.9:
            print xmlfile, str(feat_start) + ":" + str(feat_stop) + "\t" + "(" + str(aa_length), \
            "<", str(hit_len_list[0]) + ")" + "\t" + def_list[0] + "\t" + "is too short. Likely a fragment."
        elif int(aa_length) > float(hit_len_list[0]) * 1.1:
            print xmlfile, str(feat_start) + ":" + str(feat_stop) + "\t" + "(" + str(aa_length), \
            ">", str(hit_len_list[0]) + ")" + "\t" + def_list[0] + "\t" + "is too long. Check!"
        else:
            print xmlfile, str(feat_start) + ":" + str(feat_stop) + "\t" + "(" + str(aa_length), \
            "and", str(hit_len_list[0]) + ")" + "\t" + def_list[0], "are within range!"
else:
    print "No hit"

result_handle.close()
